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Abstract- In  this  paper,  we  propose  a  method  of  solving  the 
biomagnetic  inverse  problem  consisting  of  two  approaches,  the 
first  of  which  is  pointwise  normalization.  In  the  conventional 
normalization  technique,  each  variable  is  normalized  individu¬ 
ally.  This  variablewise  normalization  is  appropriate  for  scalar 
fields,  but  not  for  vector  fields  where  one  vector  on  a  grid  point 
is  represented  by  several  variables.  Hence,  pointwise  operation 
is  needed  for  vector  fields.  The  second  is  a  combination  of  norms. 
Use  of  the  i 2-norm  as  the  cost  function  of  an  optimization  prob¬ 
lem  is  known  to  lead  to  spatially  spread  solutions,  while  the 
f  i -norm  leads  to  sparse  solutions.  To  control  sparseness  of  solu¬ 
tions,  we  propose  to  use  an  internal  division  of  the  f2-norm  and 
pointwise  normalized  (/-norm.  The  optimization  problem  con¬ 
structed  as  above  can  be  recast  as  a  second-order  cone  program 
(SOCP),  a  nonlinear  convex  problem.  The  problem  can  be 
solved  using  recently  developed  efficient  interior-point  methods. 
Computer  simulations  showed  that  the  sparseness  of  estimators 
obtained  with  the  proposed  method  reflects  both  the  ratio  of 
internal  division  and  the  sparseness  of  true  sources.  Regulariza¬ 
tion  of  normalization  and  relaxation  of  constraint  conditions  in 
the  presence  of  noise  are  also  presented. 
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zation 

I.  Introduction 

Estimation  of  the  source  current  distribution  from  the  meas¬ 
ured  biomagnetic  field  is  called  a  biomagnetic  inverse  prob¬ 
lem  and  has  received  much  attention  for  many  years.  One  of 
these  problems  is  the  MEG  (magnetoencephalography)  in¬ 
verse  problem,  which  is  to  reconstruct  the  source  current  dis¬ 
tribution  from  the  magnetic  field  generated  by  the  brain. 

The  problem  is  usually  underdetermined  such  that  the  so¬ 
lution  is  not  unique.  One  approach  to  determining  the  solu¬ 
tion  uniquely  is  to  formulate  an  optimization  problem.  An 
example  of  cost  functions  of  optimization  problems  is  the 
//-norm  of  the  solution.  The  minimal  /2-norm  solution  can  be 
obtained  analytically  by  using  the  pseudoinverse  matrix,  and 
is  known  to  be  spread  spatially.  Another  example  is  the 
/[-norm.  In  this  case,  the  problem  is  solved  by  using  linear 
programming.  The  minimal  A -norm  solution  is  known  to  be 
sparse. 

Normalization  of  variables  was  proposed  for  the  use  of  the 
/rnorm  [1],  The  normalization  technique  is  operated  vari¬ 
ablewise  on  the  source  vector  and  is  appropriate  for  scalar 
fields.  There  is,  however,  room  for  taking  into  account  that 
current  distributions  are  vector  fields.  The  problem  formula¬ 
tion  and  these  conventional  methods  are  described  in  Section 
II. 

In  Section  III,  we  propose  a  new  normalization  method, 
pointwise  normalization.  This  can  be  successfully  applied  to 
vector  fields.  In  addition,  regularization  of  normalization  is 
presented.  We  also  propose  a  cost  function  that  is  an  internal 


division  of  the  //-norm  and  the  pointwise  normalized  £r norm. 
Computer  simulation  results  of  the  proposed  method  are 
shown  in  Section  IV.  In  Section  V,  a  consideration  of  noise  is 
included. 

II.  Preliminaries 

A.  System  Equation 

The  relationship  between  a  biomagnetic  field  and  a  current 
density  distribution  is  expressed  as  the  linear  equation 

m  =  Lq  ,  (1) 

where  wieRM,  q  e  Rv  ,  and  Le  RMx3Ar  are  the  meas¬ 
ured  magnetic  field  data  vector  ( M  :  number  of  sensors), 
source  current  distribution  of  discretized  points  in  the  source 
region  ( 3 N  :  number  of  points  including  orthogonal  coordi¬ 
nate  system),  and  leadfield  matrix  (or  transfer  matrix),  re¬ 
spectively.  Since  the  number  of  unknown  variables  is  greater 
than  the  number  of  equations,  M  <N  ,  the  inverse  problem 
will  be  underdetermined.  Here,  let  us  assume  rank[L]  =  M  . 
In  (1),  q  is  the  unknown  source  to  be  estimated. 

B.  Inverse  Solution 

Representing  a  solution  as  a  map  g  and  the  resulting  es¬ 
timator  as  a  vector  q  ,  we  obtain 

q=g(m)=g(Lq).  (2) 

For  example,  any  linear  inverse  solution  g  can  be  ex¬ 
pressed  by  the  matrix  G  e  R3A,xM  as 

q  =Gm  ,  (3) 

and  thus  (2)  can  be  written  explicitly. 

However  g  is  not  necessarily  written  as  a  function  with  an 
explicit  form.  An  example  is  the  case  where  the  estimator  is 
given  as  the  solution  of  an  optimization  problem. 

g(»t)=  argmin{s(<jr)|»i  =  Lq  }  (4) 

9 

s  (q)  is  the  cost  function  to  be  minimized  under  constraint 
(1). 

C.  Performance  Criterion 

One  way  to  design  a  solution  g  is  to  introduce  a  per¬ 
formance  criterion  on  g  .  Here  we  consider  the  square  sum 
of  estimation  errors  for  fundamental  vectors, 

3iV 

E(g, L,{e,})=  £||g(Le,)-e,||2 .  (5) 

1=1 

E  >  0  holds  for  any  linear  g  if  the  problem  is  underde- 
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termined.  G  (a  linear  g )  minimizing  (5)  is  the  pseudoin¬ 
verse  matrix  of  L : 

G  =  L+  =Lt{lLtY  .  (6) 

The  important  point  to  observe  is  that  the  same  estimator 
as  (6)  is  obtained  by  the  following  optimization  problem: 

g(»i)  =  argmin{  |qr|2  Int  =  Lq  }.  (7) 

Although  a  linear  g  satisfying  E  =  0  yields  the  true  so¬ 
lution,  it  cannot  occur  due  to  the  singularity  of  L  in  fl). 
Generally  speaking,  a  nonlinear  g  can  be  focused  more 
point  sources  through  the  minimization  of  E,  whereas  the 
estimators  for  their  superposition  are  usually  not  the  superpo¬ 
sition  of  the  estimators  for  each  point  source. 

D.  Normalized  i j-norm 

Matsuura  and  Okabe  [1]  proposed  in  an  optimization 
problem  employing  a  variablewise  normalized  frnorm  as  the 
cost  function: 

g[(;w)  =  D  '  argmin{s,(^)  \m  =  Lq  },  (8) 

where  L  =  LD  1 ,  q  =  Dq  ,  and  D  is  a  diagonal  matrix 
whose  diagonal  terms  are  |/1||2,...,|Z3Ar|o  .  lt  is  the  i  th 

column  vector  of  L  .  Thus  we  have  ||/(.  ||  '  and 

qi  =  \lt\\q: ,  i.e.,  the  variablewise  normalization. 

The  cost  function  S,  (r/  )  in  (8)  is 

3  N 

Si(?)=||flf||1 =LINL  •  (9) 

;=i 

The  optimization  problem  can  be  formulated  as  a  linear 
program.  Solutions  are  known  to  be  sparse  with  their  number 
of  nonzero  variables  being  equal  to  or  smaller  than  M  if 
L  is  nondegenerate.  Particularly,  for  point  sources,  g,  in 
(8)  successfully  achieves  E  =  0 . 

Note  that  the  variablewise  normalization  presented  above 
has  the  intended  affect  only  when  source  q  is  a  scalar  field. 
If  not,  normalization  that  allows  vector  fields  is  required. 
This  case  will  be  further  investigated  in  the  next  section. 


to  each  other.  This  nonorthogonality  causes  the  variablewise 
normalization  to  be  unsuitable  for  vector  fields. 

For  vector  fields,  we  propose  a  pointwise  normalization 
instead  of  the  variablewise  normalization.  That  is,  instead  of 
(8),  we  use 

g plM  =  H  l  argmin{s pl{q)\m  =  Lq  },  (11) 

q 


where  L  =  LH  1  and  q  =  Hq  .  H  is  a  block  diagonal 
matrix  given  as  follows: 

l'l  =v,Ey:' 


ii  =  E)nv? 


(12) 


H  = 


H, 


H, 


where  eigenvalue  decomposition  of  I  Ll  is  utilized.  Thus 
we  have  L,  =  TP,/!, 1/2  and  qt  =  E'nV.' q. ,  i.e.,  the  point- 


wise  normalization.  Let 


r  =  [r, , . . . ,  rN  ]  ' , 
and  the  cost  function  Spl  (q) 


(13) 


SpiWHk 


P 


LIN 

i=l 


(14) 


This  is  the  fi-norm  of  the  vector  r  ,  which  consists  of  the 
pointwise  normalized  norm  of  each  point. 

By  solving  (11)  the  following  theorem  holds. 


Theorem  Suppose  we  are  solving  (1),  or  equivalent,  (10).  Let 
L  and  L  (i  =£  j)  have  no  common  subspace  other  than  0 
(zero  vector),  i.e., 

rank[L  L;]  =  rank[L  ]  +  ranklZ^.]  =  6,  if  i  ^  j  ,  (15) 

and  each  pointwise  transfer  matrix  be  of  full  rank,  i.e., 

rank[LJ  =  3.  (16) 


Then,  any  point  source  is  correctly  estimated  by  solving  (11). 
A  point  source  with  active  point  a  is: 

[^0,  i  =  a 
1  =  0,  i  a 


<h 


(17) 


III.  Proposed  method 
A.  Pointwise  Normalization 


Denoting  the  transfer  matrix  of  point  i  by  L;  e  RM*  , 


(1)  can  be  rewritten  as 


m  =  [L, 


qN 


(10) 


where  qt  e  R3  is  the  source  vector  on  the  i  th  grid  point. 


In  the  case  of  vector  fields,  more  than  one  variables  related 
to  the  same  grid  point  altogether  represent  the  vector  on  that 
point.  Those  variables  in  general  produce  m  s  not  orthogonal 


The  proof  is  not  included  due  to  space  limitation. 

This  result  includes  E(gpl,L,{c,})=  0  and  the  fact  that 
even  if  there  are  more  than  one  nonzero  variables  in  q  ,  q 
is  correctly  solved  when  the  nonzero  variables  are  on  the 
same  discretized  point.  While  the  former  is  already  shown  in 
[1]  under  a  stronger  assumption  than  (15)  and  (16),  the  non¬ 
degeneracy  assumption,  the  latter  cannot  be  achieved  by  the 
variablewise  normalization. 

In  norm  minimizations,  a  point  that  produces  smaller 
magnetic  fields  than  other  points  tends  to  have  a  smaller  es¬ 
timated  value.  Normalization  is  intended  to  help  avoiding 
such  a  dependence  of  the  estimated  values  on  measurement 
properties.  Estimators,  however,  for  source  distributions  other 
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than  point  sources  could  be  perturbed  if  the  matrix  HJH 
has  a  large  condition  number.  We  use  a  regularization  tech¬ 
nique  to  contend  with  this  issue: 

H'  =  (D  +  a/)l/2V,T  (18) 

Spl(?,a)  =  £||//'9,.||2,  (19) 

i=i 

where  a  is  the  regularization  parameter  for  normalization. 
Estimations  are  thought  to  be  more  stable  for  larger  a , 
whereas  normalization  has  no  effect  as  a  — »  °°  . 

B.  Combination  of  Norms 

So  far,  we  have  seen  that  the  linearity  of  solution  g  and 
performance  E(g,L,{e,})  =  0  are  incompatible.  Because  g 
is  desired  to  have  both  properties,  a  compromise  between 
linearity  and  performance  for  point  sources  should  be  made. 

One  approach  is  to  linearly  combine  an  estimator  by  linear 
estimation  (6),  which  is  the  same  as  (7),  and  one  by  estima¬ 
tion  (11).  The  obtained  estimator,  however,  has  no  additional 
information  compared  to  the  two  original  estimators  and  in¬ 
deed,  seems  useless  in  practice. 

In  contrast,  using  internal  division  of  the  cost  functions  of 
(7)  and  (11)  leads  to  an  interesting  estimator.  The  new  cost 
function,  the  combined -norm,  is 

Spl2(9,a,(3)=(3||4  +  (l-(3)Spl(a,9),  (20) 

where  p,  0  <  p  <  1  is  the  ratio  of  the  internal  division.  The 
estimators  for  different  P  show  sparse  to  spread  character¬ 
istics  as  p  varies  from  0  to  1 . 

The  proposed  cost  function  (20)  is  easily  shown  to  be  a 
norm.  Thus  the  problem 

gpl2(/M,a,p)  =  argmin{spi2(9,a,p)|/M  =  L^}  (21) 

is  a  convex  optimization  problem.  Indeed,  this  problem  (21) 
is  readily  reformulated  as  a  second-order  cone  program 
(SOCP)  [2].  Although  this  reformulation  includes  the  relaxa¬ 
tion  of  limiting  conditions,  equivalence  of  the  optimization 
problem  is  preserved.  For  SOCPs,  several  efficient  inte¬ 
rior-point  methods  have  been  developed  (see  [3]  for  exam¬ 
ple.) 


Fig.  1 .  The  estimated  results  for  scalar  fields  on  a  line,  without  normalization: 
A)  a  point  source,  B)  a  Gaussian  spread  source.  1)  The  true  distribution,  2) 
the  min  ^i-norm  estimator,  3)  min  combined-norm  with  p  =  0.1,  4)  min 

combined-norm  with  p  =  0.3  ,  and  5)  min  ^-norm.  It  is  shown  that  the 
spreads  of  the  minimum  combined-norm  estimators  correspond  to  that  of  the 
true  distributions,  while  those  of  the  minimum  i\-  and  ^2-norm  estimators  do 
not. 


IV.  Results 

A.  Scalar  Fields  on  a  Line 

Figure  1  shows  estimators  for  a  point  source  and  a  spread 
source,  both  on  a  line.  The  simulation  system  is  the  same  as 
the  one  used  in  [4],  i.e.,  a  one -dimensional  model. 

In  combined-norm  cases,  the  spread  of  estimators  reflects 
the  internal  division  ratio  and  the  spread  of  the  true  distribu¬ 
tion,  while  the  spreads  of  both  the  minimum  (rnorm  and 
(2-norm  estimators  do  not. 

B.  a  Vector  Field  on  a  Plane 

Estimated  results  for  a  vector  field  on  a  plane  having  two 
active  regions  are  illustrated  in  Fig.  2.  The  plane  is  assumed 
to  be  disk  shaped  with  radius  8.06cm  and  placed  at  4cm 
height  in  a  spherical  conductor.  Magnetic  fields  are  observed 
by  64  sensors,  whose  radiuses  are  from  10.4cm  to  15.2cm. 
The  number  of  grid  points,  N  ,  is  688. 

In  the  minimization  of  the  /rnorm  and  combined-norm, 
the  regularization  of  the  pointwise  normalization  is  per¬ 
formed  such  that  the  condition  number  after  the  operation 
XmjHJH  +cd)/Xrmn{HJH  +  al)=104  ,  where  XmjA) 
and  kn.n{A)  are  the  maximum  and  the  minimum  eigenvalue 
of  matrix  A  ,  respectively.  Without  the  regularization, 
KjHTH)/Lmm{HTH)=l.06xl0\ 

For  the  distribution  shown  in  Fig.  2,  p  =  0.01  seems  appro¬ 
priate  for  approximating  the  source.  Note  that  the  effect  of 
the  mixing  parameter,  p ,  will  differ  if  the  system  transfer 
matrix  and/or  the  number  of  variables,  3 N  ,  differ. 

V.  Discussion  and  conclusions 
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Fig.  2.  The  estimated  results  for  a  vector  field  on  a  plane.  The 
darkness  of  each  dot  corresponds  to  the  norm  of  the  vector  on 
each  point.  1)  The  true  distribution,  2)  the  minimum  pointwise 
normalized  ^i-norm  estimator,  3)  min  combined-norm  with 
P  =  0.01 , 4)  min  ^2-norm.  The  normalization  is  regularized. 
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In  practice,  noise  is  inevitably  mixed  into  measurements. 
Assuming  noise  has  a  normal  distribution,  relaxation  of  con¬ 
straint  conditions  as  follows  may  be  a  suitable  approach: 

sPi2n  (?>  a,  (3,  y)  =  ySpl2  ( q ,  a,  p)  +  (l  -  y)|| m  -  Lq\\2  (22) 
Spi2n  ( w,  oc,  p,y)  —  arg  min{Spl2n  oc,  p,  y)l-  (23) 

Problem  (23)  is  still  an  SOCP.  The  second  term  on  the  right 
hand  side  of  (22)  may  be  modified  by  the  noise  covariance 
matrix. 

In  this  paper,  we  proposed  a  pointwise  normalization 
method,  which  can  solve  point  sources  correctly  when  used 
with  the  €rnorm,  to  handle  vector  fields  appropriately.  The 
regularization  of  the  normalization  is  also  introduced  in  order 
to  obtain  stable  estimations.  Then,  a  function  derived  by  in¬ 
ternal  division  of  the  G-norm  and  pointwise  normalized 
z1 1 -norm  is  proposed  as  a  cost  function  of  the  optimization 
problem.  The  problem  is  solved  by  reformulating  it  to  an 
SOCP;  the  solution  shows  sparse  or  spread  characteristics 
according  to  the  ratio  of  internal  division. 

The  next  step  in  the  research  is  to  determine  function  pa¬ 
rameters  both  in  the  case  where  one  has  some  a  priori 


knowledge  about  source  activity  and  where  nothing  is  as¬ 
sumed  previously. 
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